Create Jacobian template parameter overload for constrain functions #2559

SteveBronder · 2021-08-19T04:34:01Z

Summary

Related to stan-dev/stanc3#947 this PR allows stanc3 to pass a bool template parameter to the constraints to decide whether or not they should accumulate the jacobian calculation into the lp parameter. See the PR above to note see examples of how this is used.

All I really did here was add a new overload to the *_constrain functions that has a bool Jacobian template parameter. If that is flipped to true then we do the jacobian calculation.

Tests

I changed all of the current constrain tests to use the new Jacobian interface, which should still cover all the current test examples and the new API.

Side Effects

Release notes

Adds an overload for the constrain functions on whether to accumulate jacobians into log probability argument.

Checklist

Math issue #(issue number)
Copyright holder: Steve Bronder

The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
- Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
the basic tests are passing
- unit tests pass (to run, use: ./runTests.py test/unit)
- header checks pass, (make test-headers)
- dependencies checks pass, (make test-math-dependencies)
- docs build, (make doxygen)
- code passes the built in C++ standards checks (make cpplint)
the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested

…e supplied directly to the constrain function

…4.1 (tags/RELEASE_600/final)

stan-buildbot · 2021-08-20T03:43:07Z

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	2.9	2.96	0.98	-1.89% slower
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.98	-2.56% slower
eight_schools/eight_schools.stan	0.1	0.1	1.01	1.03% faster
gp_regr/gp_regr.stan	0.16	0.15	1.01	1.22% faster
irt_2pl/irt_2pl.stan	5.84	5.78	1.01	1.04% faster
performance.compilation	89.11	86.73	1.03	2.67% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	8.72	8.61	1.01	1.3% faster
pkpd/one_comp_mm_elim_abs.stan	30.5	30.76	0.99	-0.84% slower
sir/sir.stan	141.44	123.7	1.14	12.54% faster
gp_regr/gen_gp_data.stan	0.03	0.03	0.99	-1.45% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	3.01	2.99	1.01	0.56% faster
pkpd/sim_one_comp_mm_elim_abs.stan	0.39	0.41	0.97	-3.59% slower
arK/arK.stan	2.52	1.85	1.36	26.67% faster
arma/arma.stan	0.76	0.83	0.92	-9.06% slower
garch/garch.stan	0.67	0.71	0.94	-6.12% slower
Mean result: 1.02300453647

Jenkins Console Log
Blue Ocean
Commit hash: ca90f8c

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

bob-carpenter

I think you may fire me as reviewer after all these change requests. Most of them are around doc and I'm happy to help with that if you'd like---just let me know.

bob-carpenter · 2021-08-23T15:14:25Z

stan/math/prim/fun/cholesky_corr_constrain.hpp

+ */
+template <bool Jacobian, typename T>
+inline auto cholesky_corr_constrain(const T& y, int K, scalar_type_t<T>& lp) {
+  if (Jacobian) {


[question]
I'm unclear on how this function gets resolved since it seems to call itself inside the if (Jacobian) block. Is this only resolved because the nested call doesn't provide a template parameter?

[comment]
I opened an issue to refactor the implementation: #2561

Yes since there's no default template value for Jacobian it diverts to the signatures without the Jacobian template parameter

bob-carpenter · 2021-08-23T15:21:50Z

stan/math/prim/fun/cholesky_factor_constrain.hpp

-Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
-cholesky_factor_constrain(const T& x, int M, int N, value_type_t<T>& lp) {
+inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
+cholesky_factor_constrain(const T& x, int M, int N, scalar_type_t<T>& lp) {
  check_size_match("cholesky_factor_constrain", "x.size()", x.size(),


[optional]
Ideally the error messages will be phrased in a way that a Stan user can understand them from the Stan program. So maybe "constrained size" would be better than "x.size()". I don't think this can fail for Stan programs because we manage the sizes, and it's not part of this issue, so I'm marking "optional".

Yep just fixed these up to `constrained size'.

bob-carpenter · 2021-08-23T15:22:44Z

stan/math/prim/fun/cholesky_factor_constrain.hpp

+template <bool Jacobian, typename T>
+inline auto cholesky_factor_constrain(const T& x, int M, int N,
+                                      scalar_type_t<T>& lp) {
+  if (Jacobian) {


[optional]
I generally prefer a ternary operator when appropriate, so this would be

return Jacobian ? cholesky_factor_constrain(x, M, N, lp) : cholesky_factor_constrain(X, M, N);

if it fits or

return Jacobian ? cholesky_factor_constrain(x, M, N, lp) : cholesky_factor_constrain(X, M, N);

imo I like the if (Jacobian) for these just because it future proofs us a bit so that if we ever move to C++17 we can pretty much just do a copy replace for if (Jacobian) for if constexpr (Jacobian)

Interesting. Does the ternary operator not get reduced at compile time if the condition is constexpr? Obviously the template parameters are all constexpr!

bob-carpenter · 2021-08-23T15:26:15Z

stan/math/prim/fun/corr_constrain.hpp

+ *
+ * <p>\f$\log | \frac{d}{dx} \tanh x  | = \log (1 - \tanh^2 x)\f$.
+ *
+ * @tparam Jacobian If true, incremented `lp` with the log Jacobian


[optional]
This should just be "increment" and it's "log absolute Jacobian determinant" if you want to be precise. Or in this case, "log absolute derivative" or if you want to be verbose "log absolute derivative of the constraining transform".

I went with the bit of doc you wrote for these in the comment below

bob-carpenter · 2021-08-23T15:26:24Z

stan/math/prim/fun/corr_constrain.hpp

@@ -49,6 +49,31 @@ inline auto corr_constrain(const T_x& x, T_lp& lp) {
  return tanh_x;
 }

+/**
+ * Return the result of transforming the specified scalar or container of values


The doc for all of these needs to mention that the Jacobian is incremented with the log absolute Jacobian determinant of the transform if the template parameter Jacobian is set to true. The goal is that the first sentence of the doc should give as thorough and precise an explanation of the function's arguments and returns as possible. If you don't think it fits in one sentence, add a second sentence on the Jacobian.

Yes I appended the doc you mention below with all of these so should be good now

Looks good. Thanks.

bob-carpenter · 2021-08-23T16:15:06Z

stan/math/prim/fun/unit_vector_constrain.hpp

+
+/**
+ * Return the unit length vector corresponding to the free vector y.
+ * See https://en.wikipedia.org/wiki/N-sphere#Generating_random_points


same comments as other doc on return, etc.

refer to our reference manual, which in turn has references into the literature

bob-carpenter · 2021-08-23T16:15:39Z

test/unit/math/mix/fun/cholesky_corr_constrain_test.cpp

-  typename stan::scalar_type<T>::type lp = 0;
-  stan::math::cholesky_corr_constrain(x, inv_size(x), lp);
+auto g3(const T& x) {
+  stan::scalar_type_t<T> lp = 0;


[comment]
thanks for tidying up!

bob-carpenter · 2021-08-23T16:16:43Z

test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp

@@ -3,7 +3,8 @@
 TEST(mathMixMatFun, cholesky_factor_constrain) {
  auto f = [](int M, int N) {
    return [M, N](const auto& x1) {
-      return stan::math::cholesky_factor_constrain(x1, M, N);
+      stan::scalar_type_t<std::decay_t<decltype(x1)>> lp = 0.0;


[question]
Doesn't scalar_type_t do all the decaying for you? If not, what's this for?

Yep, fixed!

bob-carpenter · 2021-08-23T16:17:31Z

test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp

    };
  };

  auto f2 = [](int M, int N) {
    return [M, N](const auto& x1) {
      stan::scalar_type_t<std::decay_t<decltype(x1)>> lp = 0.0;
-      stan::math::cholesky_factor_constrain(x1, M, N, lp);
+      stan::math::cholesky_factor_constrain<true>(x1, M, N, lp);


[question]
Are we still testing cholesky_factor_constrain? Given that it's still in the API, it should still be tested. Especially since I've filed the issue to refactor the implementations.

Same question for all of these tests.

Yep! See my comment here

bob-carpenter · 2021-08-23T16:18:54Z

test/unit/math/mix/fun/lub_constrain_matvar_test.cpp

  };
  auto f3 = [](const auto& x, const auto& lb, const auto& ub) {
    stan::return_type_t<decltype(x), decltype(lb), decltype(ub)> lp = 0;
-    stan::math::lub_constrain(x, lb, ub, lp);
+    stan::math::lub_constrain<true>(x, lb, ub, lp);


Isn't this redundant with the test on line 12? Are we still testing the old lub_constrain? We should be.

Yes so for each of these tests corresponds to the constraint without lp and is for the constrain with lp. In order to test the API I thought it would just be easier to change all of these to use it, though eod it still differs to the original functions the test was checking

…emplate-constrains

…4.1 (tags/RELEASE_600/final)

…emplate-constrains

stan-buildbot · 2021-08-24T18:18:12Z

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	2.91	2.97	0.98	-2.18% slower
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.99	-1.17% slower
eight_schools/eight_schools.stan	0.1	0.1	0.98	-1.71% slower
gp_regr/gp_regr.stan	0.16	0.15	1.01	0.7% faster
irt_2pl/irt_2pl.stan	5.89	5.84	1.01	0.88% faster
performance.compilation	89.49	86.85	1.03	2.95% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	8.7	8.38	1.04	3.74% faster
pkpd/one_comp_mm_elim_abs.stan	30.92	28.98	1.07	6.25% faster
sir/sir.stan	125.44	120.62	1.04	3.84% faster
gp_regr/gen_gp_data.stan	0.03	0.03	1.0	-0.36% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	2.99	2.92	1.02	2.21% faster
pkpd/sim_one_comp_mm_elim_abs.stan	0.39	0.41	0.95	-4.83% slower
arK/arK.stan	1.87	1.85	1.01	1.14% faster
arma/arma.stan	0.83	0.91	0.91	-9.32% slower
garch/garch.stan	0.71	0.63	1.12	10.5% faster
Mean result: 1.01059146362

Jenkins Console Log
Blue Ocean
Commit hash: 2971378

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

bob-carpenter · 2021-08-24T22:51:17Z

Is this ready for a final review? I should be able to get to it tomorrow (Wednesday).

SteveBronder · 2021-08-24T23:01:06Z

Yep! Though I had two doc qs for the below

#2559 (comment)

…constrains

stan-buildbot · 2021-08-26T00:36:36Z

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	3.47	3.51	0.99	-1.13% slower
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	1.02	1.81% faster
eight_schools/eight_schools.stan	0.09	0.09	1.02	2.41% faster
gp_regr/gp_regr.stan	0.14	0.14	0.97	-3.45% slower
irt_2pl/irt_2pl.stan	5.11	5.17	0.99	-1.17% slower
performance.compilation	91.37	89.29	1.02	2.27% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	8.21	8.11	1.01	1.2% faster
pkpd/one_comp_mm_elim_abs.stan	31.61	30.73	1.03	2.79% faster
sir/sir.stan	125.34	119.83	1.05	4.4% faster
gp_regr/gen_gp_data.stan	0.03	0.03	0.96	-4.59% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	3.14	2.99	1.05	4.91% faster
pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.4	0.94	-5.92% slower
arK/arK.stan	2.04	2.75	0.74	-34.87% slower
arma/arma.stan	0.25	0.28	0.9	-10.77% slower
garch/garch.stan	0.77	0.68	1.14	11.97% faster
Mean result: 0.988603611096

Jenkins Console Log
Blue Ocean
Commit hash: e4650cc

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

bob-carpenter

Looks good. Thanks! I left some residual comments, but no to-do items.

bob-carpenter · 2021-08-26T19:54:31Z

stan/math/prim/fun/cholesky_factor_constrain.hpp

+template <bool Jacobian, typename T>
+inline auto cholesky_factor_constrain(const T& x, int M, int N,
+                                      scalar_type_t<T>& lp) {
+  if (Jacobian) {


Interesting. Does the ternary operator not get reduced at compile time if the condition is constexpr? Obviously the template parameters are all constexpr!

bob-carpenter · 2021-08-26T19:54:51Z

stan/math/prim/fun/corr_constrain.hpp

@@ -49,6 +49,31 @@ inline auto corr_constrain(const T_x& x, T_lp& lp) {
  return tanh_x;
 }

+/**
+ * Return the result of transforming the specified scalar or container of values


Looks good. Thanks.

bob-carpenter · 2021-08-27T16:24:54Z

stan/math/prim/fun/unit_vector_constrain.hpp

 *  scalar type.
 * @param y vector of K unrestricted variables
 * @return Unit length vector of dimension K
 */
 template <typename T, require_eigen_col_vector_t<T>* = nullptr,
          require_not_vt_autodiff<T>* = nullptr>
-inline auto unit_vector_constrain(const T& y) {
+inline plain_type_t<T> unit_vector_constrain(const T& y) {
  using std::sqrt;
  check_nonzero_size("unit_vector_constrain", "y", y);


No. The idea is that y[N] = sqrt(1 - sum(y[1:N-1])). So if N = 1, the sum is zero, and we're left with sqrt(1). That's good, because [1]' is a unit vector, but [0]' is not!

I opened the new issue here: #2569

SteveBronder and others added 2 commits August 18, 2021 16:31

Add a new overload for the constrain functions so that Jacobian can b…

dfa27aa

…e supplied directly to the constrain function

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

ca90f8c

…4.1 (tags/RELEASE_600/final)

bob-carpenter mentioned this pull request Aug 23, 2021

single templated implementation of constraints #2561

Open

4 tasks

bob-carpenter requested changes Aug 23, 2021

View reviewed changes

SteveBronder and others added 5 commits August 23, 2021 20:47

Merge remote-tracking branch 'origin/develop' into feature/jacobian-t…

d74af94

…emplate-constrains

Update docs for constrain functions with Jacobian template parameter

661ce32

fix corr doc

e9797e1

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

318dd13

…4.1 (tags/RELEASE_600/final)

Merge remote-tracking branch 'origin/develop' into feature/jacobian-t…

2971378

…emplate-constrains

Merge remote-tracking branch 'origin' into feature/jacobian-template-…

e4650cc

…constrains

bob-carpenter approved these changes Aug 27, 2021

View reviewed changes

bob-carpenter merged commit 9532191 into develop Aug 27, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Create Jacobian template parameter overload for constrain functions #2559

Create Jacobian template parameter overload for constrain functions #2559

SteveBronder commented Aug 19, 2021

stan-buildbot commented Aug 20, 2021

bob-carpenter left a comment

bob-carpenter Aug 23, 2021

SteveBronder Aug 23, 2021

bob-carpenter Aug 23, 2021

SteveBronder Aug 24, 2021

bob-carpenter Aug 23, 2021

SteveBronder Aug 24, 2021

bob-carpenter Aug 26, 2021

bob-carpenter Aug 23, 2021

SteveBronder Aug 24, 2021

bob-carpenter Aug 23, 2021

SteveBronder Aug 24, 2021

bob-carpenter Aug 26, 2021

bob-carpenter Aug 23, 2021

bob-carpenter Aug 23, 2021

bob-carpenter Aug 23, 2021

SteveBronder Aug 24, 2021

bob-carpenter Aug 23, 2021

SteveBronder Aug 24, 2021

bob-carpenter Aug 23, 2021

SteveBronder Aug 24, 2021

stan-buildbot commented Aug 24, 2021

bob-carpenter commented Aug 24, 2021

SteveBronder commented Aug 24, 2021

stan-buildbot commented Aug 26, 2021

bob-carpenter left a comment

bob-carpenter Aug 26, 2021

bob-carpenter Aug 26, 2021

bob-carpenter Aug 27, 2021

Create Jacobian template parameter overload for constrain functions #2559

Create Jacobian template parameter overload for constrain functions #2559

Conversation

SteveBronder commented Aug 19, 2021

Summary

Tests

Side Effects

Release notes

Checklist

stan-buildbot commented Aug 20, 2021

bob-carpenter left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stan-buildbot commented Aug 24, 2021

bob-carpenter commented Aug 24, 2021

SteveBronder commented Aug 24, 2021

stan-buildbot commented Aug 26, 2021

bob-carpenter left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment